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Abstract 
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1  Project  objectives 

The  specific  research  objectives  of  this  research  were  proposed  as: 

1.  To  develop  novel  computational  statistical  inference  approaches  for  dynamic  vehicle  state  es¬ 
timation,  through  a  combination  of  machine  learning  and  reduced-order  modeling  techniques. 

2.  To  develop  adaptive  model  reduction  approaches  that  use  dynamic  data  to  update  reduced- 
order  models  for  vehicle  flight  limit  prediction. 

3.  To  develop  approaches  for  online  management  of  multifidelity  models  and  sensor  data,  using 
variance-based  sensitivity  analysis. 

4.  To  apply  our  new  methods  to  quantify  the  benefits  of  a  self-aware  UAV  in  terms  of  reliability, 
maneuverability  and  survivability. 

The  project  has  accomplished  all  these  objectives.  Objective  3  was  achieved  using  support 
vector  machines  and  self-organizing  maps,  rather  than  variance-based  sensitivity  analysis.  The 
remainder  of  this  report  summarizes  the  main  project  outcomes  and  accomplishments. 

2  Project  Outcomes  and  Accomplishments 

2.1  Summary 

The  project  has  resulted  in  the  development  of  new  DDDAS  methodology  and  DDDAS  algorithms, 
new  models  for  a  DDDAS-enabled  self-aware  unmanned  aerial  vehicle  (UAV),  and  a  demonstration 
of  the  value  of  DDDAS  in  the  context  of  dynamic  data-driven  structural  assessment  to  support 
decision-making  for  a  damaged  vehicle  taking  evasive  action  in  a  hostile  environment.  The  specific 
project  accomplishments  include  the  following: 

•  A  conceptual  design  framework  for  a  DDDAS-enabled  UAV.  The  framework  includes 
sizing  of  the  UAV  at  the  overall  vehicle  level  and  detailed  finite  element  analysis  at  the  panel 
level.  It  also  includes  a  structural  response  model  that  incorporates  multiple  degradation  or 
failure  modes  including  damaged  panel  strength  (BVID,  thru- hole),  damaged  panel  stiffness 
(BVID,  thru-hole),  loose  fastener,  fretted  fastener  hole,  and  disbonded  surface. 

•  A  new  data-driven  approach  for  the  online  updating  of  the  flight  envelope  of  a 
DDDAS-enabled  self-aware  UAV  subjected  to  structural  degradation.  The  main 
contribution  of  this  part  of  the  work  is  a  general  methodology  that  leverages  both  physics- 
based  modeling  and  data  to  decompose  tasks  into  two  phases:  expensive  offline  simulations 
to  build  an  efficient  characterization  of  the  problem,  and  rapid  data-driven  classification  to 
support  online  decision-making.  In  the  approach,  physics-based  models  at  the  wing  and  vehi¬ 
cle  level  run  offline  to  generate  libraries  of  information  covering  a  range  of  damage  scenarios. 
These  libraries  are  queried  online  to  estimate  vehicle  capability  states.  The  state  estima¬ 
tion  and  associated  quantification  of  uncertainty  are  achieved  by  Bayesian  classification  using 
sensed  strain  data. 

•  A  demonstration  of  the  value  of  DDDAS  on  a  conceptual  UAV  executing  a  pull-up 
maneuver,  where  the  vehicle  flight  envelope  is  updated  dynamically  with  onboard 
sensor  information.  During  vehicle  operation,  the  maximum  maneuvering  load  factor  is 


estimated  using  structural  strain  sensor  measurements  combined  with  physics-based  informa¬ 
tion  from  pre-computed  damage  scenarios  that  consider  structural  weakness.  Compared  to 
a  baseline  case  that  uses  a  static  as-designed  flight  envelope,  the  self-aware  vehicle  achieves 
both  an  increase  in  probability  of  executing  a  successful  maneuver  and  an  increase  in  overall 
usage  of  the  vehicle  capability. 

•  A  new  data-driven  strategy  for  online  structural  assessment  at  the  panel  level  for 
a  DDDAS-enabled  self-aware  UAV.  During  the  offline  phase,  high-fidelity  finite  element 
simulations  are  used  to  construct  reduced-order  models  and  classification  criteria:  proper 
orthogonal  decomposition  approximations  and  self-organizing  maps  are  combined  to  realize  a 
fast  mapping  from  measured  quantities  to  system  capabilities.  During  the  online  phase,  the 
surrogate  mapping  is  employed  to  directly  estimate  the  vehicle’s  evolving  structural  capability 
from  sensor  data.  The  approach  has  been  demonstrated  for  a  test  problem  of  a  composite  wing 
panel  on  an  unmanned  aerial  vehicle  that  undergoes  degradation  in  structural  properties. 

•  A  new  algorithm  for  fast  kernel  density  estimation  in  high  dimensions.  We  intro¬ 
duced  novel  methods  for  pruning  and  approximating  the  far  field.  Our  far  field  approximation 
requires  only  kernel  evaluations  and  does  not  use  analytic  expansions.  Pruning  is  not  done 
using  bounding  boxes  but  rather  combinatorially  using  a  sparsified  nearest-neighbor  graph  of 
the  input.  The  time  complexity  of  our  algorithm  depends  linearly  on  the  ambient  dimension. 
The  error  in  the  algorithm  depends  on  the  low-rank  approxinrability  of  the  far  field,  which  in 
turn  depends  on  the  kernel  function  and  on  the  intrinsic  dimensionality  of  the  distribution  of 
the  points.  The  error  of  the  far  field  approximation  does  not  depend  on  the  ambient  dimen¬ 
sion.  We  report  results  for  Gaussian  kernel  sums  for  100  million  points  in  64  dimensions,  for 
one  million  points  in  1000  dimensions,  and  for  problems  in  which  the  Gaussian  kernel  has  a 
variable  bandwidth.  To  the  best  of  our  knowledge,  all  of  these  experiments  are  impossible  or 
prohibitively  expensive  with  existing  fast  kernel  summation  methods. 

•  Fast  algorithms  for  Bayesian  inverse  medium  problems.  We  consider  a  supervised 
inverse  medium  problem  algorithm  using  a  Bayesian  framework  and  a  variational  formulation 
for  a  maximum  a  posteriori  (MAP)  estimation  of  the  label  field.  In  the  Bayesian  framework, 
we  must  define  the  likelihood  function  (the  probability  of  the  data  given  the  label  field)  and 
the  prior  probability  (for  the  label  field).  We  propose  a  non-parametric,  high-dimensional, 
kernel  density  estimation  (KDE)  for  the  likelihood  function,  based  on  Gabor  features  (300 
per  pixel).  This  approach  better  approximates  non-local  correlations.  For  the  prior  function, 
we  use  a  simple  smoothness  prior.  We  provide  experimental  evidence  that  this  likelihood 
function  performs  very  well  and  converges  to  the  correct  segmentation  with  the  number  of 
training  datasets. 

These  accomplishments  are  described  in  more  detail  in  the  following  sections. 

2.2  Aircraft  Design 

A  conceptual  design  of  a  UAV  was  established  to  allow  for  the  assessment  of  structural  damage 
and  degradation  cases  on  realistic  maneuvers  from  a  vehicle.  The  design  was  completed  using 
first-principles  sizing  and  Federal  Aviation  Regulation  (FAR)  23  guidelines  (similar  to  NASA  CR- 
2010-216794/VOL2  and  AFOSR  Final  Report  for  contract  FA9550-10-C-0175).  Aurora  utilized 
an  internally  developed  Bayesian  Multi-Disciplinary  Optimization  (BMDO)  code  to  generate  the 
conceptual  design.  The  vehicle  has  a  wing  span  of  55  feet,  a  cruise  speed  of  140  knots  at  an  altitude 


of  25,000  feet,  and  a  500  pound  payload  capability.  The  range  of  the  aircraft  is  estimated  to  be 
roughly  2500  nautical  miles,  corresponding  to  a  duration  of  17.5  hours.  An  illustration  of  the 
vehicle  is  shown  in  Figure  1.  The  estimated  weights  of  the  component  comprising  the  conceptual 
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Figure  1:  Conceptual  aircraft  design. 


design  are  included  in  Table  1  and  the  relative  percentage  of  each  component  is  shown  in  Figure  2. 
Characteristic  speeds  of  the  conceptual  UAV  are  included  in  Table  2  and  characteristics  of  the 
lifting  surfaces  of  the  vehicle  are  included  in  Table  3.  The  design  of  this  conceptual  UAV  allows  for 
adequate  range  and  duration  to  explore  capability  as  a  function  of  the  changing  structural  state  of 
the  vehicle. 

The  aerodynamic  characteristics  of  the  wings  are  determined  using  the  low  Reynolds  number 
airfoil  design  and  analysis  code,  XFOIL.  XFOIL  is  an  interactive  program  for  the  design  and 
analysis  of  subsonic  isolated  airfoils.  The  Reynolds  number,  Re,  is  a  dimensionless  number  that 
gives  a  measure  of  the  ratio  of  inertial  forces  to  viscous  forces: 


Re  = 


vL 


(1) 


v 

where  v  is  the  mean  velocity  of  the  aircraft,  L  is  a  characteristic  linear  dimension,  and  v  is  the 
kinematic  viscosity.  Reynolds  numbers  are  use  in  airfoil  design  to  manage  “scale  effect”  when 
computing/comparing  characteristics.  For  this  work,  the  mean  velocity  is  taken  as  the  flight  speed, 
V,  the  characteristic  length  is  the  average  chord  length,  c,  and  the  kinematic  viscosity  of  the  air  is 
taken  at  cruise  altitude: 


Re 


Re  = 


Vc 


2.126  x  10~5  — 

S 

(72.0^)  (1.614m)  , 

- - ^ =  5.47  x  106 

2.126  X  10“5?2U 
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Table  1:  Component  weights  of  the  conceptual  UAV 


Component 

Weight  (lbs.) 

Fuselage 

476.2 

Wing 

286.6 

Payload 

500.0 

Vertical  Tail 

22.4 

Horz.  Tail 

7.5 

Engine 

423.5 

Pylon 

25.0 

Fuel 

1060.8 

HPE  System 

225.0 

Nose  LG 

31.4 

Main  LG 

174.1 

Sum 

3232.5 

Main  Landing  Gear  Weight  (5%) 


Engine  Weight  (13%)  Pylon  Weight  (1%) 


Figure  2:  Pie  chart  of  the  estimated  weights  of  the  components  comprising  the  conceptual  UAV 
design. 


Results  from  the  XFOIL  analysis  are  displayed  in  Figures  3  and  4,  the  lift  coefficient  versus  the 
angle  of  attack  and  drag  coefficient,  respectively.  In  order  to  approximate  the  lift  at  cruise,  the 
wing  lift  coefficient,  Cl,  is  assumed  equal  to  the  airfoil  lift  coefficient,  Cf. 

W  =  L  =  qSCL  =  qSCt  (4) 


where  the  airfoil  coefficient  is  calculated  for  140  KTAS  at  15000  feet  to  be: 


Q 


1 

q 


1  /3232.5 lbs\ 

41.6^  V  146.7 /f2  ) 


0.5291 


(5) 


The  resulting  wing  incidence  required  for  level  fuselage  during  cruise  is  0.885°.  The  plots  of  Cl 
and  Cl/Cd  versus  the  angle  of  attack  are  shown  in  Figure  5.  The  maximum  lift  coefficient  (clean) 
depends  on  the  airfoil  characteristics.  The  following  equations  are  valid  for  high-aspect-ratio  wings 
with  moderate  sweep: 

Lmax  =  qSCLrnax  (6) 


Table  2:  Characteristic  speeds  of  the  conceptual  UAV 


Aircraft  Speed 

Speed  [KTAS] 

Speed  [m/s] 

VA 

66.9 

34.4 

VS 

57.5 

29.6 

Vc 

(cruise  and  loiter) 

140.0 

72.0 

VD 

168.0 

86.4 

Table  3:  Characteristic  of  the  conceptual  UAV 


Metric 

Variable 

Value 

Takeoff  Weight 

W 

3257.3 

lbs 

Wing  Area 

S 

146.7 

in2 

Wing  Sweep 

A 

4.3 

O 

Wing  Span 

b 

54.2 

ft 

Wing  Root  Chord 

Co 

3.79 

ft 

Wing  Tip  Chord 

Ct 

1.52 

ft 

Wing  Aspect  Ratio 

AR 

20 

Front  Spar  Location 

ifspar 

15 

% 

Rear  Spar  Location 

^ir  spar 

70 

% 

Wing  Airfoil 

DA01 

Horz.  Tail  Airfoil 

NACA  0012 

Vert.  Tail  Airfoil 

NAVA  0012 

CLmax  =0.9Cimax  cos  Ac  =  0.9  (1.8263)  cos  (4.3°)  =  1.64 


Lr 


=  ^41.6  ^  (146.7  ft2)  (1.64)  =  10002  lbs 
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10002  lbs 


^ m.ax  — 


=  3.09 
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W  3232.5  lbs  ^ 

This  load  factor  is  less  than  3.8,  which  would  correspond  to  a  CLmax  of  2.01.  Note,  these  calcu¬ 
lations  do  not  account  for  control  surfaces.  A  plot  from  Raymer’s  Aircraft.  Design:  A  Conceptual 
Approach  [11]  provides  maximum-lift  trends  versus  sweep  angle  for  several  classes  of  aircraft,  shown 
in  Figure  6.  From  the  wing  loading  and  maximum  lift  coefficient,  the  stall  speed  of  the  aircraft  can 
be  determined: 

W  =  L  =  qstallSCLmax  (10) 


(11) 


Vstall  —  \  /  I  n 


where  the  clean  wing  stall  speed  at  various  altitudes  is: 


Qstall  —  1 

Therefore,  the  stall  velocity  is: 


W\  1 


S)  cLn 


f 3232.5  lbs 
V  146.7  ft2 

I  ^qstall 


1  ,lbs 
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(12) 


/ 1  Q'l 


Figure  3:  Lift  coefficient  versus  angle  of  attack  for  Re  =  5.47  x  106. 


Figure  4:  Lift  coefficient  versus  drag  coefficient  for  Re  =  5.47  x  106. 


The  stall  velocity  at  three  different  altitudes  of  interest  can  be  calculated: 

Vstali  (h  =  0 ft)  =  106.3—  =  63.0  knots  (14) 

s 

Vstau  (h  =  15000 ft)  =  134.3-  =  79.6  knots  (15) 

s 

Vstaii  (h  =  25000 ft)  =  159.3-  =  94.4  knots  (16) 

s 

The  stall  speed  at  15000  feet,  Equation  15,  is  below  the  designed  140  knots  (72.0  m/s)  cruise  speed, 
as  indicated  in  Table  2.  This  means  that  the  wing  will  not  stall  at  cruise,  but  the  stall  speed  is  less 
than  the  desired  57.5  knots  stall  speed.  The  lack  of  control  surfaces  increases  stall  speed: 


Vmi  ■  v  l  ( s )  czL  k  {cL,  (l7) 

Mandatory  separation  from  stall  speed,  as  specified  in  FAR  23,  pushes  the  minimum  maneuver 
speed  above  the  design  point: 

vAmin  =  VsVn 


(18) 
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Figure  5:  Lift  coefficient  and  lift  coefficient  divided  by  drag  coefficient  versus  angle  of  attack. 
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Table  4  contains  the  minimum  speeds  for  the  three  altitudes. 


(20) 


Table  4:  Minimum  maneuver  speeds,  as  mandated  in  FAR  23,  for  the  conceptual  design 


Aircraft  Speed 

Speed 

Altitude  =  0  ft 

Altitude  =  15000  ft 

Altitude  =  25000  ft 

KEAS 

KTAS 

m/s 

KEAS 

KTAS 

m/s 

KEAS 

KTAS 

m/s 

vArni„ 

122.8 

122.8 

63.2 

155.1 

195.8 

100.8 

183.9 

275.6 

141.8 

Vs„t„ 

63.0 

63.0 

32.4 

79.6 

100.5 

51.7 

94.4 

141.4 

72.7 

VCmi„ 

154.4 

154.4 

79.4 

155.1 

195.8 

100.8 

184.0 

275.6 

141.8 

VDmin 

215.9 

215.9 

111.1 

216.9 

273.9 

140.9 

257.3 

385.4 

198.3 

For  initial  assessment  of  the  DDDAS  algorithms,  the  aircraft  performs  a  standard  rate  turn.  A 
sustained  standard  rate  turn  is  performed  by  banking  the  aircraft  at  an  angle,  a,  and  increasing 
the  lift  normal  to  the  wing  planform.  Increasing  the  lift  is  necessary  to  balance  the  aircraft  weight 
while  creating  a  net  centripetal  force  perpendicular  to  the  velocity,  V,  of  the  aircraft.  The  constant 
turn  radius,  r,  is  equal  to  the  velocity,  V,  divided  by  the  turn  rate, 


r  = 


V 

ip 


(21) 


The  turn  rate,  ip,  is  calculated  as: 


ip  = 


gV n2  —  1 

V 


(22) 


where  g  is  the  acceleration  of  gravity,  n  is  the  load  factor  (ratio  of  lift  of  the  aircraft  to  its  weight), 
and  V  is  the  velocity.  The  variables  involved  with  the  standard  turn  are  illustrated  in  Figure  7. 
Equation  22  can  be  rearranged  to  give  the  load  factor  as  a  function  of  aircraft  velocity  and  turn 
rate:  _ 


n  = 


N 


(23) 


Performing  a  standard  rate  turn  allows  the  aircraft  to  change  its  heading  at  a  constant  velocity;  a 
larger  bank  angle  allows  for  a  higher  rotational  speed  and  a  tighter  radius  of  turn,  while  increasing 
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Wings  of  moderate  aspect  ratio  (4-8) 


Source:  Raymer  Aircraft  Design:  A  Conceptual  Approach 

Figure  6:  Maximum  lift  coefficient  versus  angle  of  attack  [11], 


the  lift  force  necessary  to  keep  the  aircraft  at  altitude.  Increasing  the  load  factor  increases  the 
loading  on  the  structure  of  the  aircraft. 

The  lift  distribution  along  the  wing  for  multiple  turn  rates  and  multiple  radii  is  estimated 
using  Schrenk’s  approximation  based  on  an  elliptical  and  chord-wise  estimation  of  lift.  This  lift 
distribution  is  reacted  by  the  wing  box  structure  of  the  concept  UAV  in  bending  and  shear.  The 
cross-section  of  an  example  wing  box  is  shown  in  Figure  8.  The  wing  spar  caps  react  the  majority  of 
bending  while  the  wing  spar  webs  react  the  majority  of  shear.  These  load  reactions  induce  tension 
and  compression  stress  in  the  lower  and  upper  wing  skins,  respectively,  due  to  bending;  shear  stress 
is  induced  in  the  spar  web  due  to  the  shear  loading.  The  basic  relations  for  calculating  the  shear 
stress  in  the  web  and  the  tensile  or  compressive  stress  in  the  cap  are  given  in  Equations  24  and  25. 


'Tweb  — 


5 -*-(>?) 
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■^web(rf)  C±  (?))  twebs  (hf  spar  T  hr  spar') 
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Figure  7:  Standard  rate  turn  variables. 


Figure  8:  Cross-section  of  an  example  wing  box. 
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These  relations  allow  for  estimating  the  turn  rate  or  radius  of  the  concept  UAV  based  on  the 
capability  of  the  structural  components  to  react  the  loading.  As  shown  in  Figure  8,  a  constant- 
percent  chord  control  surfaces  can  be  placed  in  front  of  or  behind  the  wing  box.  Trailing  edge  flaps 
can  occupy  the  last  30%  of  the  chord  and  increase  the  lift  coefficient  and  increase  the  wing  area 
(Fowler-type  extending  flaps).  The  increase  in  lift  coefficient  is: 

C 

A  Ci  =  1.9—  (26) 

Leading  edge  slats  can  occupy  the  first  15%  of  the  chord  and  prevent  premature  airflow  separation 
caused  by  the  flap. 

r" 

A  Ci  =  0.4—  (27) 

Further  investigating  the  increase  in  the  coefficient  of  lift  due  to  the  addition  of  flaps,  the  change 
in  the  lift  coefficient  can  be  calculated  via  Equation  28: 

A CLmax  =  0.9A Clmax  (  SI^A  cosAh.l.  (28) 

\  bref  J 

Estimation  of  the  maximum  increase  in  lift  coefficient,  assume  the  control  surfaces  run  the  entire 
span,  as  shown  in  Figure  9  of  the  wing  and  that  the  leading  edge  extends  the  chord  by  5%: 


A  CLmax  =  0. 9(0.4)  (1.05)cos(6.4°)  =  0.38 


and  assume  the  trailing  edge  extends  the  chord  by  10%: 


A  CL 


max 


(29) 


0.9(1.9)(l.l)cos(3.1°)  =  1.87 


(30) 


The  total  maximum  lift  coefficient  becomes  3.89  (1.64  +  0.38  +  1.87). 

The  maximum  lift  and  stall  speed  should  now  be  revisited.  Recalculating  the  maximum  lift  and 
load  factor: 


L 


max 


^41.6  ^  (146.7  ft2)  (3.89) 

_  Lmax  _  23740  lbs  _ 
W  3232.5  lbs 


=  237401bs 

7.34 
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(32) 


Recalculating  the  stall  speeds: 


Qstall  — 


W\  1 

s'/ 


f 3232.5  lbs 
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1  lbs 

- =  5.66—^ 

3.89  ft2 


Vgtaii  (h  =  0 ft)  =  69.0—  =  40.9 knots 

s 

ft 

Vstau  (h  =  15000 ft)  =  87.2—  =  51.6 knots 

s 

ft 

Vstau  (h  =  25000 ft)  =  103.4—  =  61.3 knots 

s 


(33) 

(34) 

(35) 
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From  the  updated  stall  speeds,  V-n  diagrams  are  created  for  the  three  altitudes.  These  V-n  diagrams 
are  shown  in  Figures  10  through  12.  The  updated  minimum  speeds  at  the  three  altitudes  are 


Flight  Envelope  at  Altitude=0ft 


Figure  10:  Flight  envelope  at  an  altitude  of  0  feet. 

updated  in  Table  5.  Vhmrri  is  less  than  the  designed  140  knots  cruise  speed.  Therefore,  the  aircraft 
is  able  to  maneuver  at  speeds  of  140  knots  (KTAS)  at  15000  feet  (minimum  maneuver  speed  is 
127.1  knots).  While  the  cruise  speed  of  the  initially  designed  UAV  was  too  low  to  meet  FAR  23 
due  to  wing  loading,  the  necessary  lift  coefficient  at  15000  feet  was  reduced  by  50%.  The  viable 
design  is  therefore  used  in  moving  forward  with  the  example  problem. 


Flight  Envelope  at  Altitude=  15 000ft 


Figure  11:  Flight  envelope  at  an  altitude  of  15000  feet. 


Table  5:  Minimum  maneuver  speeds,  as  mandated  in  FAR  23,  for  the  conceptual  design  with  flaps 


Aircraft  Speed 

Speed 

Altitude  =  0  ft 

Altitude  =  15000  ft 

Altitude  =  25000  ft 

KEAS 

KTAS 

m/s 

KEAS 

KTAS 

m/s 

KEAS 

KTAS 

m/s 

79.7 

79.7 

41.0 

100.7 

127.1 

65.4 

119.4 

178.8 

92.0 

vSmt„ 

40.9 

40.9 

21.0 

51.6 

65.2 

33.5 

61.3 

91.8 

47.2 

vCmi„ 

154.4 

154.4 

79.4 

154.4 

194.9 

100.3 

154.4 

231.3 

119.0 

VnmiT, 

215.9 

215.9 

111.1 

215.9 

272.7 

140.3 

215.9 

323.5 

164.4 

2.3  Finite  Element  Modeling 

The  magnitude  of  wing  loading  on  the  concept  UAV  is  directly  related  to  the  turn  rate;  therefore 
a  section  of  the  wing  is  used  as  an  example  panel  on  which  to  test  the  algorithms  developed  under 
the  DDDAS  program.  This  allows  for  damage  or  degradation  that  reduces  the  structural  capability 
of  the  panel  to  react  a  sustained  load  without  failure  to  directly  affect  the  ability  of  the  concept 
UAV  to  perform  a  standard  rate  turn. 

In  order  to  reduce  the  magnitude  of  loading  in  anticipation  for  experimental  testing,  a  section  of 
the  wing  is  chosen  outboard  at  a  location  where  the  wing  skin  is  sized  for  strength  (further  outboard 
the  wing  skin  is  sized  by  minimum  gage  or  handling  loads  rather  than  flight  loading) .  At  a  location 
roughly  260  inches  outboard  along  the  wing,  the  wing  skin  consisting  of  4  plies  of  MTM45-1/AS4 
carbon  composite  in  a  quasi-isotropic  layup  is  necessary  to  meet  the  loading  requirements  due  to 
FAR  23  design  requirements,  while  not  exceeding  strain  design  allowable  values  for  the  plies.  For 
the  finite  element  models,  a  “cut-out”  section  of  the  wing  is  modeled.  The  relative  location  of  this 
“cut-out”  along  the  span  of  the  wing  is  shown  in  Figure  13. 

Three  structural  conditions  were  initially  considered  for  the  example  panel:  pristine,  moderate, 
and  severe  damage.  Damage  cases  are  defined  as  delamination  of  the  first  ply,  which  is  typical  of  a 
low-velocity  impact  on  the  panel  or  interlaminar  separation  caused  by  stress  concentrations  under 


Flight  Envelope  at  Altitude=25000ft 


Figure  12:  Flight  envelope  at  an  altitude  of  25000  feet. 


cyclic  loading  due  to  embedded  foreign  object  debris  or  other  defects.  The  panel  conditions  are 
modeled  using  the  finite  element  method  and  strain  values  for  the  panel  are  estimated  for  a  variety 
of  loading  conditions  that  can  be  traced  back  to  loading  from  a  variable  turn  rate.  The  level  of 
delamination  and  resulting  strain  fields  for  the  three  conditions  are  shown  in  Figure  14.  The  images 
of  the  panel  finite  element  model  show  the  extent  and  location  of  the  varying  delamination  size, 
where  the  red  and  light  orange  color  indicates  delamination  between  the  first  and  second  plies.  The 
panels  were  reinforced  around  the  boarders  of  the  panel  with  through  holes  to  simulate  mounting 
the  wing  skins  to  the  structure.  The  panel  predicted  strains  indicate  strain  gradients  that  are 
expected  to  be  present  within  the  panel  for  the  given  loading.  These  sets  of  strains  can  then  be 
used  by  the  vehicle  system  to  identify  the  condition  of  the  panel  based  on  measured  strain  values 
from  strain  gauges  or  other  sensors  on-board  the  vehicle.  Once  the  panel  condition  is  determined, 
the  capability  of  the  vehicle  to  perform  a  sustained  turn  necessary  to  follow  a  flight  path  can  be 
determined  and  a  go/no-go  decision  can  be  made. 

The  material  properties  used  within  the  finite  element  models  are  provided  by  the  material 
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Figure  13:  Relative  location  of  the  “cut-out”  panel  along  the  span  of  the  wing. 
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Figure  14:  Finite  element  models  and  resulting  strain  fields  for  the  three  structural  conditions. 


manufacturer  (www.umeco.com).  To  account  for  uncertainty  within  the  finite  element  models,  a 
normal  distribution  and  coefficient  of  variance  (Cv)  of  5%  is  assumed  and  values  corresponding  to 
99.86%  confidence  (3cr)  are  used.  A  summary  of  the  properties  is  provided  here: 

•  Composite  properties: 

—  Advantages:  high  strength/weight  ratio;  directional  strength  /  stiffness  tailoring;  fatigue 
resistant 

—  Disadvantages:  more  complex;  brittle;  most  design  allowables  /  testing  is  proprietary 
—  MTM45-1/CF0525-36%RW,  193gsm  3K  AS4  plain  weave  carbon  fiber 

*  Matrix:  MTM45-1  flexible  curing  temperature,  high  performance,  toughened  epoxy 
matrix  system  optimized  for  low  pressure,  vacuum  bag  processing 

*  Fiber:  AS4  carbon  fiber  filaments  made  from  PAN  (polyacrylonitrile)  precursor, 
surface  treated;  3,000  filaments/tow 

*  36%  Resin  Weight  (54.34%  Fiber  Volume) 

*  193  grams  per  square  meter  (gsm)  areal  weight  (0.201  mm  ply  thickness) 

*  Plain  Weave  warp  and  weft  criss-cross  in  over-under  pattern 

2.4  Structural  Response  Model 

A  structural  response  model  was  created  for  the  conceptual  UAV  using  ASWING.  ASWING  is 
a  program  for  the  aerodynamic,  structural,  and  control-response  analysis  of  aircraft  with  flexible 
wings  and  fuselages  of  high  to  moderate  aspect  ratio.  The  program  couples  the  structural  and 
aerodynamical  response  of  the  vehicle.  Thus,  the  loading  on  the  aircraft  wing  can  be  investigated  for 
different  flight  conditions.  MATLAB  scripts  link  the  finite  element  models  discussed  in  Section  2.3 
with  the  ASWING  models.  The  degradation  of  mechanical  properties  results  in  aerodynamic 
changes  (i.e. ,  shape  of  loaded  wing  will  be  different  due  to  changes  in  structural  response).  The 


coupled  models  allow  for  vehicle  response  to  be  investigated.  The  results  from  these  coupled  models 
indicate  that  the  capability  of  the  UAY  will  change  between  2  and  12%  (corresponding  to  strain 
levels  within  the  panel  changing  between  1  and  20%)  for  the  levels  of  structural  damage  investigated. 

The  model  includes  multiple  degradation  or  failure  modes  including:  damaged  panel  strength 
(BVID,  thru- hole),  damaged  panel  stiffness  (BVID,  thru- hole),  loose  fastener,  fretted  fastener  hole, 
and  disbonded  surface.  The  algorithm  is  implemented  in  object-oriented  MATLAB  coding  using 
handle  classes  for  speed  and  compatibility  with  the  following  algorithmic  elements: 

•  Framework  (main. nr) 

—  Maintains  “truth”  data  regarding  flight,  environment,  and  events  (inputs) 

•  Inputs  (Environment. m,  MeasuredDanrage.m,  VehicleManeuver.m) 

—  Measures  inputs  using  sensors  or  flight  controller;  should  include  signal  uncertainty, 
noise,  and  fault /failure  models 

—  Can  be  replaced  by  hardware-in-the-loop 

•  Algorithm  (VehicleCapability.m) 

—  Interrogates  input  classes  and  operates  on  the  information  to  provide: 

=i=  diagnostic  information  (Monitor  method); 

*  current  or  future  vehicle  state  information  (State  method); 

*  airframe  capability  (Capability  method);  and 

*  prognostic  information  (Eval  method) 

Structural  Modeling  Routine 

Tasks  2-4  are  iterative. 

1.  Initialize 

•  Read  finite  element  model  (nas  file) 

•  Generate  nominal  node  map  and  element  strength  capability 

2.  Monitor  (damage, environment) 

•  Receive  measured  damage  information  (type,  size,  location) 

•  Receive  temperature  information  (gross  temperature) 

•  Modify  local  strength  capability 

•  Modify  local  stiffness  or  remesh  (not  implemented  yet) 

3.  State(nraneuver) 

•  Export  and  run  updated  finite  element  model 

•  Import  finite  element  model  stresses 

•  Set  current  maneuver  state 


4.  Capability 


•  Calculate  finite  element  model  failure  indices  against  load  cases 

•  Notch  load  cases  based  on  reduced  capability  (not  implemented  yet) 

5.  Evaluate 

•  Estimate  Remaining  Useful  Life  (RUL)  using  simulated  load  spectrum  (not  implemented 
yet) 

•  Cost  function  information  to  decision  making  algorithms  (not  implemented  yet) 

Multiple  diagnostic  models  can  be  integrated  into  the  Monitor  method.  Error  checking  mod¬ 
els  can  be  integrated  into  the  State  method  to  compare  updated  model  state  predictions  and 
additional  sensors.  Multiple  prognostic  models  can  be  integrated  into  the  Eval  method,  e.g., 
Sendeckyj /Whitworth  residual  strength  models  for  composite  fatigue  and  composite  damage  pro¬ 
gression  models  for  overloading  conditions.  A  load  spectrum  generated  using  multiple  loading 
scenarios  generated  by  decision  making  agent  to  determine  cost  of  short-  and  long-term  operation. 
Current  state  and  R.UL  estimates  are  then  combined  with  updated  usage  scenarios. 

2.5  Experimental  Testing 

Experimental  testing  was  planned  to  follow  model  development.  This  was  an  additional  activity 
that  was  not  part  of  the  original  proposal;  however,  it  was  added  because  if  the  potential  value 
in  validating  our  modeling  and  algorithm  contributions.  As  of  the  writing  of  this  report,  multiple 
experimental  specimens  containing  different  types  of  damage  have  been  manufactured  and  a  custom 
made  load  frame  is  under  development.  Testing  of  the  specimens  will  occur  as  follow-on  work  at  a 
later  time. 

The  mechanical  aspects  of  the  load  frame  are  complete,  with  all  components  assembled.  Addi¬ 
tional  software  development  is  required  in  order  to  control  the  load  frame.  The  original  goal  was 
to  have  a  table-top  experiment.  In  order  for  the  frame  to  supply  loads  large  enough  to  damage  the 
composite  panels,  the  frame  evolved  to  be  a  free-standing,  self-reacting,  load  frame,  as  shown  in 
Figure  15.  The  frame  is  still  of  a  size  that  it  can  be  transported  and  assembled  on-site  (unlike  typ¬ 
ical  mechanical  testers  that  require  permanent  installation).  Quasi-dynamic  tensile  loading  allows 
the  loading  experienced  by  the  panel  at  the  specified  location  on  the  wing  and  for  specific  vehicle 
maneuvers  to  be  simulated.  The  load  frame  can  apply  uniaxial  tensile  loading.  The  “cut-out” 
section  of  the  wing  was  selected  such  that  uniaxial  tension  is  a  good  approximation  of  the  actual 
loading  experienced.  The  magnitude  of  loading  is  to  be  directly  related  to  executed  maneuver.  The 
software  drivers  will  match  the  applied  loading  to  anticipated  loading  resulting  from  maneuvers  of 
the  aircraft. 

The  load  frame  use  a  Duff-Norton  worm  gear  actuator  to  impart  load  into  the  system.  The 
actuator  is  rated  for  5  tons  (10000  lbf)  and  is  driven  via  a  DC  motor.  The  motor  is  a  Leeson  21 
amp,  12  V-DC  drive.  An  Elmo  motion  controller  (DC- WHISTLE  20/60)  controls  the  DC  motor 
and  interfaces  with  a  personal  computer  in  order  to  send  control  signals.  A  Lab  VIEW  VI  is  under 
development  to  provide  the  motor  control  signals  as  well  as  record  the  data  from  the  load  cell  and 
displacement  sensor.  The  control  will  permit  either  load  control  or  displacement  control  based 
on  real-time  data  from  the  load  cell  or  the  LVDT,  respectively.  Controller  feedback  is  provided 
via  encoder  (attached  to  actuator  worm  gear),  LVDT,  and/or  load  cell.  The  position  sensor  is  a 
MacroSensors  DC750-500  LVDT.  The  load  cell  is  a  Sensor  Development  10188-014,  10  kip  threaded 
rod  load  cell  (uniaxial  tensile  load).  The  encoder  is  an  Accu-Coder.  The  flight  profile  determined 
from  the  vehicle  path  planner  will  set  the  load  path  (load  history).  The  load  path  will  be  executed 
in  real-time,  with  real-time  sensor  data  provided  to  the  multi-fidelity  models. 


Figure  15:  Load  frame  with  composite  specimen  and  DAQ  system. 


A  National  Instruments  cDAQ  system  with  multiple  modules  provides  input/output  control 
between  the  computer  and  the  test  instrumentation.  One  aspect  that  will  be  investigated  before 
experimental  testing  is  started  is  integrating  optical  strain  measurements  in  order  to  measure  full 
field  strains,  as  opposed  to  single  point  measurements  via  strain  gages. 

The  test  specimens  were  manufactured  at  Aurora’s  Mississippi  plant.  Layups  of  four  large  panels 
were  fabricated  and  sectioned  into  a  total  of  200  test  specimens.  While  manufacturing  these  panels, 
two  types  of  damage  were  included.  Controlled  areas  of  delamination  were  created  by  including 
teflon  shapes  between  the  first  and  second  plies,  as  shown  in  Figure  16.  The  second  form  of  damage 
was  simulated  fiber  fracture  that  was  implemented  by  cutting  lengths  of  fibers  within  a  single  ply. 
The  specimens  were  autoclave  cured.  An  image  of  the  panels  placed  within  the  autoclave  after  the 
cure  cycle  is  shown  in  Figure  17. 


Figure  16:  MIT  graduate  student  Marc  Lecerf  places  teflon  inserts  between  plies  to  create  delam¬ 
ination. 


Figure  17:  Panels  being  cured  in  Aurora’s  autoclave. 


2.6  Multifidelity  DDDAS  Methods 

Figure  18  shows  our  approach  that  combines  offline  computation  with  online  sensor  data  to  provide 
a  time-constrained,  updated  estimate  of  UAV  flight  capability.  More  specifically,  we  combine 
information  from  physics-based  models,  simulated  offline  to  build  a  scenario  library,  together  with 
dynamic  sensor  data  in  order  to  estimate  current  flight  capability.  We  also  use  dynamic  data  to 
modify  our  information-gathering  strategies  to  manage  uncertainty  in  our  capability  estimates. 


Figure  18:  To  achieve  a  self-aware  vehicle  capability,  we  use  offline  physics-based  modeling  com¬ 
bined  with  dynamic  sensor  data  to  achieve  dynamically  updated  estimates  of  vehicle  capabilities 
(e.g.,  the  vehicle  flight  envelope). 


2.6.1  Problem  setup 

We  consider  two  kinds  of  quantities  of  interest: 

•  The  quantities  that  are  measured  during  flight,  referred  to  as  the  measured  quantities  of 
interest.  We  denote  these  by  qm(x),  m  =  1, . . . ,  M,  where  qm  is  the  mth  measured  quantity 
of  interest  (generally  a  vector  representing  a  discretized  field  quantity)  and  we  consider  M 
such  quantities. 

•  The  quantities  employed  in  the  decision  process  that  give  information  about  capability  and 
performance  constraints,  referred  to  as  the  capability  quantities  of  interest.  We  denote  these 
by  sc(x),  c  =  1, . . . . ,  C,  where  sc  is  the  cth  capability  quantity  of  interest  (generally  a  vector 
representing  a  discretized  field  quantity)  and  we  consider  C  such  quantities. 

Both  quantities  of  interest  are  functions  of  the  vehicle  state,  x,  which  includes  quantities  de¬ 
scribing  the  current  damage  state  (e.g.,  damage  location,  extent  and  depth).  In  an  offline  phase,  we 
use  high-fidelity  models  to  obtain  detailed  information  about  the  system  and  its  possible  responses 
for  different  flight  conditions  and  different  damage  states.  The  simulated  information  is  stored  in  a 
damage  library  and  is  used  to  build  a  mapping  from  measured  quantities  of  interest  to  capability 
quantities  of  interest  via  surrogate  models.  In  the  online  phase,  these  surrogate  models  provide 
rapid  estimates  of  vehicle  capability  given  real-time  measurements.  The  next  two  subsections  de¬ 
scribe  in  more  detail  the  models  and  methods  employed  in  our  offline  and  online  phases. 

2.6.2  Offline  stage 

Our  offline  phase  considers  a  suite  of  multifidelity  physics-based  models,  at  varying  levels  of  physics 
resolution,  and  ranging  from  the  panel  level  to  the  full  vehicle  level.  We  store  offline  analysis 
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Figure  19:  The  representation  of  our  concept  UAV  within  ASWING.  The  structure  is  specified  as 
a  set  of  interconnected  slender  beams,  where  lifting  surfaces  have  additional  aerodynamic  properties 
specified  along  their  span.  The  plots  to  the  right  show  the  aeroelastic  trim  solution  computed  by 
ASWING  for  a  pull-up  maneuver  at  a  fixed  flight  velocity  and  load  factor. 


information  in  a  damage  library.  We  use  a  variety  of  strategies  to  build  surrogate  models  from 
this  information,  including  support  vector  machines,  proper  orthogonal  decomposition  (POD)  and 
self-organizing  maps. 

Vehicle-level  models.  At  the  vehicle  level,  a  combination  of  ASWING  [3]  and  the  Variational 
Asymptotic  Beam  Cross-Sectional  Analysis  (VABS)  [2]  provides  estimates  of  internal  wing  stresses 
and  deflections  as  a  function  of  input  aircraft  kinematic  states  and  estimates  of  damage  to  the 
nominal  aircraft  structure. 

ASWING  [3]  is  a  nonlinear  aero-structural  solver  for  flexible-body  aircraft  configurations  of 
high  to  moderate  aspect  ratio.  ASWING  uses  unsteady  lifting  line  aerodynamics  together  with 
a  finite  difference  Euler-Bernoulli  beam  structural  model.  The  aerodynamic  model  is  extended 
with  Prandtl-Glauert  compressibility  treatment  and  a  sectional  stall  model.  The  nonlinear  Euler- 
Bernoulli  beam  model  permits  analysis  of  large  deflections.  Figure  19  shows  the  ASWING  represen¬ 
tation  of  our  concept  UAV  and  the  corresponding  aerostructural  assessment  for  an  example  pull-up 
maneuver.  The  ASWING  model  is  a  set  of  interconnected  slender  beams — one  each  for  the  wing, 
fuselage,  horizontal  stabilizer,  and  vertical  stabilizer.  Lifting  surfaces  (the  wing  and  stabilizers) 
have  additional  cross-sectional  lifting  properties  that  are  pre-specified. 

VABS  [2]  models  the  wing  box  as  an  array  of  two-dimensional  cross-sectional  finite  element 
models.  The  cross-sections  capture  the  details  of  a  multi-ply  composite  wing  skin  and  local  damage 
effects.  VABS  computes  lumped  stiffness  and  inertial  properties  at  a  reference  point  in  each  cross 
section,  forming  a  global  line  representation  of  the  beam.  A  standard  beam  problem  solver  finds 
the  global  force  and  moment  distribution  along  this  reference  line  given  input  forces  and  moments. 
Using  the  reference  line  solution,  the  internal  strain  field  can  be  recovered  in  the  beam  cross  sections 
using  relations  initially  computed  by  VABS.  Figure  20  provides  an  overview  of  the  VABS  modeling 
framework.  In  this  work,  we  use  the  specific  UM/VABS  implementation  developed  in  FORTRAN 


UM/VABS  calculates 

a)  Stiffness  properties  along  reference  line 

b)  Influence  coefficients  relating  reference 


Specify  2D  FEMs  line  solution  to  cross-sectional  warping 

of  cross  sections  Internal  stress 


Any  1 D  beam  solver  can  find  the 
force  +  moment  solution  along 
reference  line  based  on  boundary 
conditions 


Figure  20:  Variational  Asymptotic  Beam  Cross-Sectional  Analysis  (VABS)  allows  for  dimensional 
reduction  of  an  expensive  three-dimensional  beam  solution  into  two-dimensional  finite  element 
models  coupled  with  an  external  beam  solver. 


Figure  21:  Panel  layout  and  layer  sequence.  Panel  state  variables  include  the  damage  location, 
damage  size,  and  load  definition. 


by  R.  Palacios  and  C.  Cesnik  at  the  University  of  Michigan  [10]. 

In  our  integrated  modeling  setup,  ASWING  manages  the  one-dimensional  beam  solution,  com¬ 
puting  loads  in  the  wing  box  for  specific  flight  conditions,  while  VABS  resolves  local  stiffness  loss 
due  to  damage  on  the  aircraft  wing. 

Panel-level  models.  At  the  panel  level,  we  use  a  finite  element  model  that  simulates  and  ana¬ 
lyzes  the  panel  behavior  under  specified  loading  and  damage  scenarios  using  a  plate  model.  The 
composite  panel  comprises  carbon  fiber  plies  with  a  specified  stacking  sequence.  Four  clamped 
edges  define  boundary  conditions  that  simulate  the  presence  of  fastening  bolts  along  the  panel 
perimeter.  Figure  21  shows  a  panel  layout  and  corresponding  panel  cross-section  for  a  case  with 
four  plies  with  the  symmetric  stacking  sequence  [±45°,  0°/90°,  0°/90°,  ±45°].  The  border  area  that 
hosts  the  holes  is  reinforced  with  two  additional  plies  with  orientation  0°/90°  and  ±45°.  The  grey 
shaded  area  on  the  right  plot  of  Figure  21  denotes  the  damaged  region,  which  is  centered  at  (yd,  Zd ) 
and  has  extent  Ay  x  Az. 


The  loads  l  on  the  panel  are  defined  by  the  prescribed  aircraft  maneuver.  The  presence  of 
damage  is  simulated  by  weakening  the  stiffness  properties  of  the  finite  elements  that  belong  to 
the  prescribed  damaged  area.  This  model  is  used  to  explore  different  damage  scenarios.  For  each 
scenario,  we  generate  a  snapshot  set  by  evaluating  the  measured  quantities  of  interest  and  the 
capability  quantities  of  interest.  The  measurement  quantities  of  interest  are  the  components  of 
strain  evaluated  for  each  element  in  the  computational  mesh.  The  capability  quantities  of  interest 
are  the  failure  indices,  also  evaluated  for  each  element  in  the  computational  mesh.  Failure  index 
is  an  indicator  of  the  structural  condition  that  is  translated  into  a  scaling  factor  for  maneuver 
parameters.  It  is  defined  as  the  ratio  between  the  experienced  stress  and  the  maximum  allowable 
stress  (typically  the  compression/tension/shear  strength  that  characterizes  the  material  properties). 

Damage  library.  The  purpose  of  the  damage  library  is  to  store  information  from  analyses  run 
offline  using  the  panel- level  and  vehicle-level  models.  This  information  may  be  accessed  by  online 
analyses  in  a  variety  of  different  ways  (e.g.,  directly  accessing  the  predictions  associated  with 
a  damage  scenario,  or  using  library  entries  to  perform  dynamic  adaptation  of  a  reduced-order 
model).  The  damage  library  is  populated  by  offline  analysis  of  a  range  of  damage  scenarios  and 
kinematic  states  relevant  to  the  planned  mission  of  the  vehicle.  Information  regarding  the  planned 
mission  can  be  used  to  determine  the  type  of  damage  that  may  occur  and  the  types  of  maneuvers 
that  may  be  performed.  This  information,  along  with  resource  constraints  on  the  library,  can 
be  used  to  optimize  the  contents  of  the  damage  library  given  mission  level  objectives.  In  our 
current  implementation,  the  library  contains  model  predictions  of  sensed  strain  and  maximum 
failure  indices  for  each  scenario  considered.  Ongoing  work  is  investigating  the  utility  of  including 
model  predictions  of  other  quantities. 

Classification  and  surrogate  modeling.  We  embody  the  information  generated  by  these 
physics-based  models  using  various  surrogate  modeling  techniques.  At  the  panel  level,  we  build  low 
dimensional  representations  of  the  measured  and  capability  quantities  of  interest,  using  the  POD. 
We  then  use  self-organizing  maps  to  cluster  the  resulting  data  and  to  build  a  surrogate  model  that 
maps  measurements  to  capability  estimates.  The  steps  in  this  process  are  outlined  in  Figure  22. 
The  details  of  the  panel-level  approach  are  described  in  [8] . 

At  the  vehicle  level,  we  use  the  maximum  computed  failure  index  to  classify  all  simulated 
conditions  into  safe  (maximum  failure  index  <  1)  and  unsafe  (maximum  failure  index  >  1).  We 
build  a  representation  of  the  failure  boundary  (maximum  failure  index  =  1)  as  a  function  of  damage 
parameters  and  vehicle  operating  conditions.  Our  current  implementation  uses  a  simple  bisection 
method  to  define  the  failure  boundary;  future  work  will  use  adaptive  sampling  with  support  vector 
machines  as  in  [4].  The  steps  in  this  process  are  outlined  in  Figure  23  and  are  described  in  more 
detail  in  [7]. 

2.6.3  Online  stage 

In  the  online  stage,  our  goal  is  to  move  from  sensed  data  to  updated  estimates  of  structure  capability 
at  the  local  panel  level  and  updated  flight  envelopes  at  the  vehicle  level.  The  key  is  to  achieve  these 
estimates  sufficiently  rapidly  to  support  dynamic  decision-making,  while  leveraging  the  rich  amount 
of  physics-based  information  contained  within  our  damage  library.  We  also  aim  to  characterize  our 
level  of  confidence  in  our  updated  estimates. 

At  the  panel  level,  the  low-dimensional  POD  representations  combined  with  the  self-organizing 
map  clustering  achieve  the  desired  rapid  assessment.  Figure  22  shows  the  online  flow  of  analysis 
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Figure  22:  Diagram  of  the  offline  (black  boxes)  and  online  (red  boxes)  steps  of  our  approach  at 
the  panel  level.  POD  models  and  clusters  are  obtained  offline  and  employed  online  as  indicated  by 
the  gray  arrows. 
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Figure  23:  Diagram  of  the  offline  (black  boxes)  and  online  (red  boxes)  steps  of  our  approach 
at  the  vehicle  level.  Damage  library  information  is  stored  using  offline  clustering  based  on  the 
capability  quantities  of  interest.  The  library  information  is  employed  online  using  a  Bayes  classifier 
as  indicated  by  the  gray  arrow. 


that  moves  from  sensed  strain  data  to  estimates  of  panel  failure  indices.  The  reconstruction  of  the 
POD  coefficients  in  the  second  step  is  achieved  using  the  gappy  POD  approach  [1,5]. 

At  the  vehicle  level,  we  use  a  Bayes  classifier  that  assigns  a  posterior  estimate  of  the  probability 
of  being  in  any  damage  state  contained  in  our  damage  library.  A  simple  approach,  described  in  [7], 
assigns  equal  prior  probability  to  each  damage  state  in  the  library  and  classifies  the  vehicle’s  damage 
state  as  the  state  in  the  library  with  the  maximum  posterior  probability.  A  more  sophisticated 
approach  incorporates  mission  information,  current  maneuver,  and  a  mixture  model  of  the  posterior 
damage  state.  The  additional  mission  and  maneuver  information  permit  us  to  relax  the  assumption 
of  equal  prior  damage  state  probabilities.  The  mixture  model  of  the  posterior  damage  state  enables 
interpolation  among  records  in  the  damage  library,  permitting  us  to  characterize  a  richer  range  of 
vehicle  damage  states. 

2.6.4  Results 

The  offline  vehicle  model  allows  for  determination  of  the  vehicle  capability  over  a  range  of  damage 
events  using  classification.  Figure  24  shows  an  example  of  this  process.  We  simulate  a  pull-up 
maneuver  at  constant  velocity  V  =  260  ft/s  with  varying  load  factor  n  €  [1,  3.5]  over  a  range  of  two 
chordwise  damage  size  parameters,  with  fixed  spanwise  size  parameters  and  constant  depth.  Each 
damage  and  maneuver  combination  is  classified  as  “safe”  or  “unsafe”  based  on  the  value  of  the 
maximum  failure  index  in  the  wing  structure.  Using  this  classifier  variable,  the  boundary  surface 
between  safe  and  unsafe  maneuver  regions  for  any  given  damage  case  can  be  determined  using  a 
bisection  algorithm. 

Figure  25  presents  results  for  an  example  case  at  the  panel  level.  The  surrogate  models  are 
obtained  offline  using  an  evaluation  set  of  3000  different  damage  cases  for  fixed  panel  layout  and 
loading  condition  (Figure  21).  We  show  here  results  for  the  case  of  a  3.5220  x  7.7540  square-inch 
damage  located  at  (yd  =  10.5900,  Zd  =  12.9100)  and  involving  plies  4,  5  and  6.  The  figure  shows 
the  original  finite  element  solution  for  the  failure  index  field  over  the  panel.  We  use  this  solution  to 
generate  synthetic  data  with  which  to  test  our  approach.  The  bottom  right  plot  of  Figure  25  shows 
the  corresponding  online  estimate  of  the  failure  index  given  by  our  surrogate  modeling  approach, 
reconstructed  using  strain  sensor  measurements  that  cover  50%  of  ply  4.  The  other  two  plots  in 
Figure  25  are  included  to  give  insight  to  the  errors  due  to  the  POD  approximation  and  the  self¬ 
organizing  map  as  follows.  The  top  right  plot  represents  the  best  approximation  of  the  original 
finite  element  solution  in  the  POD  basis — this  would  be  achieved  only  if  we  could  reconstruct  the 
POD  modal  coefficients  exactly.  The  bottom  left  plot  shows  that  a  small  amount  of  additional  error 
is  introduced  by  using  the  self-organizing  map  to  map  from  POD  measurement  coefficients  to  POD 
capability  coefficients.  The  final  plot  in  the  lower  right  also  includes  the  error  due  to  inferring  the 
POD  measurement  coefficients  from  sparse  data.  Overall,  the  reconstructions  of  the  failure  indices 
are  sufficient  to  provide  a  first-cut  assessment  for  online  decision-making. 

We  also  present  representative  results  at  the  vehicle  level.  One  potential  application  of  our 
method  is  for  missions  in  contested  environments,  where  threats  to  the  vehicle  due  to  hostile  agents 
require  a  fast,  defensive  reaction  to  avoid  dangerous  regions  of  the  flight  zone.  In  addition,  the 
vehicle  may  sustain  damage  on  the  wing  surface  that  impedes  its  ability  to  operate  at  its  initial 
design  capability.  Figure  26  presents  a  schematic  of  this  scenario,  where  the  vehicle  initiates  the 
evasive  action  at  an  airspeed  of  210  ft/s  and  an  initial  load  factor  of  1.3  (representative  of  a  nominal 
maneuvering  speed  and  an  upper  bound  on  the  nominal  maneuvering  load  factor  during  normal 
operation  while  navigating  a  sequence  of  waypoints). 

We  demonstrate  how  the  decision  strategy,  informed  by  the  capability  estimate,  trades  off 
between  survivability  and  full  utilization  of  the  vehicle  capability.  For  details  on  these  results, 
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Figure  24:  Surface  representing  the  vehicle  capability  as  a  function  of  two  damage  parameters 
using  the  offline  vehicle  model.  A  pull-up  maneuver  with  constant  velocity  V  and  varying  load  factor 
n  produces  wing  loading  cases  that  are  classified  as  safe  or  unsafe,  producing  a  maximum  safe  load 
factor  Umax  for  each  damage  case  using  a  bisection  algorithm.  Note  that  when  the  chordwise  damage 
extent  is  zero,  the  damage  region  is  degenerate  with  zero  volume,  and  the  model  operates  using  the 
baseline  “pristine”  aircraft  configuration;  hence,  we  see  no  change  to  the  maximum  allowable  load 
factor. 
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Figure  25:  Failure  index  FI  over  the  panel.  The  panel  border  containing  the  bolt  holes  is  excluded 
for  plotting  clarity.  Top  left:  the  “truth”  finite  element  model  (FEM)  solution.  Top  right:  the 
representation  of  the  FEM  solution  in  the  POD  basis.  Bottom  left:  POD  approximation  combined 
with  self-organizing  map  (SOM)  surrogate,  mapping  measurement  POD  coefficients  to  capability 
POD  coefficients.  Bottom  right:  our  full  approach  using  gappy  POD  to  infer  POD  coefficients 
from  measured  data,  and  SOM  to  map  from  measurement  POD  coefficients  to  capability  POD 
coefficients. 


see  [6,7].  We  quantify  survivability  by  evaluating  probability  of  maneuver  success  over  all  possible 
damage  cases.  We  quantify  utilization  using  the  average  ratio  between  the  vehicle’s  operational 
load  factor  and  the  vehicle’s  true  maximum  load  factor. 

The  probability  of  maneuver  success,  p  (MS),  is  defined  as  the  probability  that  the  agent  chooses 
a  value  of  operating  load,  nop,  that  is  less  than  the  maximum  vehicle  load  factor  The  utiliza¬ 

tion  is  defined  as  the  average  ratio  between  the  vehicle’s  operational  load  factor  and  the  vehicle’s 
true  maximum  load  factor.  Figure  27  plots  the  probability  of  maneuver  success,  p  (MS),  versus 
the  utilization,  nutil,  for  the  two  decision  strategies:  a  static  estimate  based  on  the  vehicle  design 
flight  envelope,  and  our  DDDAS  strategy  that  uses  a  dynamic  estimate  of  the  vehicle  structural 
capability  based  on  dynamic  strain  data.  A  third  curve  is  plotted  that  shows  the  performance  of  a 
hypothetical  estimator  that  knows  the  damage  case  with  absolute  certainty;  that  is,  its  only  error 
is  due  to  the  surrogate  model  approximation  of  the  corresponding  capability  set  (in  this  case  using 
a  probabilistic  support  vector  machine). 

The  ideal  decision  strategy  would  have  both  perfect  usage  of  available  capability  (i.e.,  nutil  =  1) 
and  certain  maneuver  success  (i.e.,  p  (MS)  =  1);  this  is  marked  as  the  “Utopia”  point  in  the  upper 
right  corner  of  Figure  27.  A  sample  on  the  plot  is  considered  to  be  non-dominated  if  no  other 
sample  has  both  a  higher  p  (MS)  and  higher  n  til  value  (or  higher  value  of  one  and  equal  value  of 
the  other);  the  non-dominated  combinations  of  (nutil,  p  (MS))  for  each  estimator  are  connected  by 
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Figure  26:  Schematic  of  the  flight  scenario  to  which  we  apply  our  DDDAS  capability  estimation 
framework. 


a  dashed  line  of  the  corresponding  color. 

Figure  27  shows  that  the  static  capability  case  has  p(MS)  =  1  at  values  of  nutil  <  0.75.  This 
is  because  if  the  agent  sets  a  sufficiently  low  static  load  factor,  the  realized  load  is  less  than  . 
Thus,  within  the  scope  of  our  analysis  the  vehicle  never  exceeds  the  flight  envelope,  although  the 
loads  are  limited  to  conservative  values  and  utilization  is  low.  The  figure  also  shows  that  the  static 
capability  case  has  a  long  trail  of  samples  near  p  (MS)  =  0  at  high  values  of  nutil .  This  is  because 
at  values  of  n^p110  close  to  2.9,  the  vehicle  almost  certainly  exceeds  the  flight  envelope  unless  it  is 
in  the  pristine  case,  which  has  a  small  probability  (<  0.01)  of  occurring.  We  see  that  the  dynamic 
estimator  results  in  an  even  spread  of  points  across  the  non-dominated  front,  with  a  sharp  “knee” 
at  nutil  «  0.95  where  the  probability  of  maneuver  success  drops  rapidly. 

We  can  use  the  non-dominated  fronts  Figure  27  as  a  measure  of  performance  of  each  capability 
estimate  when  used  for  decision-making  in  the  flight  scenario.  For  instance,  if  the  agent  wants 
to  utilize  95%  of  the  maximum  vehicle  load  factor  on  average,  there  would  be  an  80%  chance 
of  maneuver  success  using  the  dynamic  estimate  of  the  load  factor  as  opposed  to  a  40%  chance 
of  success  when  operating  at  a  static  load  factor.  On  the  other  hand,  if  the  agent  can  accept 
operating  at  less  than  80%  of  the  maximum  capability  on  average,  then  both  estimators  show 
similar  performance.  We  note  this  is  most  likely  because  the  damage  cases  in  the  library  cause  a 
limited  reduction  in  the  vehicle  capability,  and  simulating  more  severe  damage  cases  would  continue 
to  emphasize  the  performance  gain  from  using  the  dynamic  capability  estimate.  Lastly,  Figure  27 
shows  that  if  the  damage  were  known  perfectly,  the  error  introduced  by  the  surrogate  modeling 
approximations  would  be  relatively  small  for  this  example.  The  difference  between  the  dynamic 
capability  curve  and  the  known  damage  curve  is  an  indication  of  the  value  (in  terms  of  vehicle 
capability  utilization  and  maneuver  success  probability)  of  increased  quality  and  quantity  of  sensors. 
Comparisons  of  this  kind  could  be  used  to  support  design  decisions  regarding  the  cost  versus  value 
tradeoffs  of  including  more  onboard  sensors. 


Figure  27:  The  probability  of  maneuver  success  (p  (MS))  versus  the  average  fraction  of  the 
vehicle  capability  utilized  {nutil)  for  an  onboard  decision  process  performed  by  an  agent  using 
either  a  static  or  dynamic  capability  estimation  strategy.  Also  plotted  is  a  hypothetical  capability 
estimate  that  knows  the  damage  with  certainty  but  still  uses  a  PSVM  capability  approximation. 
The  non-dominated  points  for  each  capability  estimation  strategy  are  connected  by  a  dashed  line 
of  the  corresponding  color. 


2.7  Fast  parallel  algorithms  for  data  assimilation  using  Bayesian  approaches 

Aircraft  state  characterization  technologies  require  fast  algorithms  for  the  solution  of  inverse  medium 
problems  for  damage  characterization.  Assuming  a  known  constitutive  law  for  the  mechanical  re¬ 
sponse  of  a  “healthy”  structure,  one  can  solve  an  inverse  medium  problem  to  find  perturbations  in 
the  material  properties  that  appear  in  the  the  constitutive  law. 

Such  problems  can  be  formulated  as  inverse  medium  problems.  We  consider  their  general  formu¬ 
lation  in  a  Bayesian  setting.  In  particular  we  considered  an  idealized  inverse  medium  problem  that 
models  the  reconstruction  of  the  unknown  structure  state  of  an  aircraft  given  partial  observations 
of  the  state.  For  that  we  considered  fast  algorithms  for  deterministic  inverse  problems,  construction 
of  likelihood  and  prior  probability  functions  for  spatial  field,  and  kernel  density  estimation. 

Below  we  discuss  some  representative  results  from  work  partially  funded  by  this  award. 

•  We  developed  a  new  algorithm  for  fast  kernel  density  estimation  in  high  dimensions  summa¬ 
rized  in  [9].  A  direct  evaluation  of  the  sum  scales  quadratically  with  the  number  of  points. 
Fast  kernel  summation  methods  can  reduce  this  cost  to  linear  complexity,  but  the  constants 
involved  do  not  scale  well  with  the  dimensionality  of  the  dataset.  The  main  algorithmic  com¬ 
ponents  of  fast  kernel  summation  algorithms  are  the  separation  of  the  kernel  sum  between 
near  and  far  field  (which  is  the  basis  for  pruning)  and  the  efficient  and  accurate  approxima¬ 
tion  of  the  far  field.  We  introduced  novel  methods  for  pruning  and  approximating  the  far 
field.  Our  far  field  approximation  requires  only  kernel  evaluations  and  does  not  use  analytic 
expansions.  Pruning  is  not  done  using  bounding  boxes  but  rather  combinatorially  using  a 
sparsified  nearest-neighbor  graph  of  the  input.  The  time  complexity  of  our  algorithm  depends 
linearly  on  the  ambient  dimension.  The  error  in  the  algorithm  depends  on  the  low-rank  ap- 
proximability  of  the  far  field,  which  in  turn  depends  on  the  kernel  function  and  on  the  intrinsic 
dimensionality  of  the  distribution  of  the  points.  The  error  of  the  far  field  approximation  does 
not  depend  on  the  ambient  dimension.  We  presented  the  new  algorithm  along  with  experi¬ 
mental  results  that  demonstrate  its  performance.  We  report  results  for  Gaussian  kernel  sums 
for  100  million  points  in  64  dimensions,  for  one  million  points  in  1000  dimensions,  and  for 
problems  in  which  the  Gaussian  kernel  has  a  variable  bandwidth.  To  the  best  of  our  knowl¬ 
edge,  all  of  these  experiments  are  impossible  or  prohibitively  expensive  with  existing  fast 
kernel  summation  methods. 

•  We  developed  fast  algorithms  for  Bayesian  inverse  medium  problems  [12].  We  consider  a 
supervised  inverse  medium  problem  algorithm  using  a  Bayesian  framework  and  a  variational 
formulation  for  a  maximum  a  posteriori  (MAP)  estimation  of  the  label  field.  In  the  Bayesian 
framework,  we  must  define  the  likelihood  function  (the  probability  of  the  data  given  the  label 
field)  and  the  prior  probability  (for  the  label  field).  In  this  work,  we  focus  on  the  likeli¬ 
hood.  Typically,  the  likelihood  of  the  intensity  is  decoupled  for  every  inverse  medium  point 
and  is  given  by  a  parametric  density  (often  a  Gaussian  function).  Instead,  we  propose  a 
non-parametric,  high-dimensional,  kernel  density  estimation  (KDE)  for  the  likelihood  func¬ 
tion,  based  on  Gabor  features  (300  per  pixel).  This  approach  better  approximates  non-local 
correlations.  For  the  prior  function,  we  use  a  simple  smoothness  prior.  We  provide  experi¬ 
mental  evidence  that  this  likelihood  function  performs  very  well  and  converges  to  the  correct 
segmentation  with  the  number  of  training  datasets. 
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Transitions 


•  Willcox  hosted  a  DDDAS  Collaboration  Meeting  at  MIT  in  May  2014.  The  meeting  brought 
together  a  small  group  of  DDDAS  researchers  with  local  industry  and  government  partici¬ 
pants,  with  a  focus  on  applications  in  the  aerospace  domain.  The  goals  of  the  meeting  were: 
(1)  to  identify  potential  beneficial  opportunities  to  transition  existing  DDDAS  methods  and 
algorithms,  including  potential  future  collaborations,  and  (2)  to  discuss  open  challenges  in 
managing  the  data-to-decisions  flow  in  future  aerospace  systems.  Government  and  industry 
participants  included  AFRL,  Lincoln  Laboratories,  Draper  Laboratories,  Raytheon,  and  Au¬ 
rora  Flight  Sciences.  The  meeting  resulted  in  productive  discussions  and  identified  several 
avenues  for  future  collaborations  among  academic  and  industry/government  participants. 

•  The  aircraft  design  framework  developed  in  part  under  this  research  has  been  used  by  de¬ 
signers  at  Aurora  Flight  Sciences  to  conduct  a  large  number  of  design  studies,  in  which  rapid 
turnaround  of  a  feasible  configuration  and  performance  estimates  were  required.  Aurora  have 
further  built  upon  these  toolsets  to  include  several  levels  of  fidelity  and  to  allow  implementa¬ 
tion  of  a  large  number  of  multidisciplinary  design  optimization  (MDO)  algorithms. 

•  Aurora  Flight  Sciences  have  an  ongoing  Air  Force  SBIR  effort  that  has  just  reached  the 
end  of  Phase  I  (we  are  working  on  Phase  II  proposal)  that  continues  to  develop/realize  the 
self-aware  vehicle  concept.  The  effort  has  focused  on  the  framework  that  connects/informs 
an  in-situ  mission  planner  of  the  capability  of  multiple  subsystems  of  a  vehicle.  Because  of 
the  DDDAS  effort,  our  technical  effort  in  this  SBIR  has  focused  on  developing  the  engine 
health/degradation  models  and  the  interface  control  with  the  mission  planner.  The  research 
has  been  directed  such  that  we  can  plug  in  the  DDDAS  structural  results  down  the  line.  The 
DDDAS  effort  is  a  key  component  necessary  for  our  framework  to  work. 


6  Supported  personnel 

Aurora  Flight  Sciences:  Jeffrey  Chambers,  David  Kordonowy. 

Massachusetts  Institute  of  Technology:  Graduate  Students:  Marc  Lecerf  (now  at  Raytheon).  Post¬ 
doctoral  Researchers:  Laura  Mainini,  Benjamin  Peherstorfer,  Demet  Ulker.  Research  Scientist: 
Douglas  Allaire  (now  at  Texas  A&M  University).  Undergraduate  Researcher:  Harriet  Li. 

University  of  Texas:  Graduate  Students:  Dhairya  Malhotra,  Bo  Xiao,  Amir  Gholami.  Postdoctoral 
Researchers:  Hari  Sundar  (now  at  University  of  Utah),  Bill  March. 
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